Optimal control of quantum superpositions in a bosonic Josephson junction 
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We show how to optimaUy control the creation of quantum superpositions in a bosonic Josephson 
junction within the two-site Bose-Hubbard model framework. Both geometric and purely numerical 
optimal control approaches are used, the former providing a generalization of the proposal of Micheli 
et al [Phys. Rev. A 67, 013607 (2003)]. While this method is shown not to lead to significant 
improvements in terms of time of formation and fidelity of the superposition, a numerical optimal 
control approach appears more promising, as it allows to create an almost perfect superposition, 
within a time short compared to other existing protocols. We analyze the robustness of the optimal 
solution against atom number variations. Finally, we discuss to which extent these optimal solutions 
could be implemented with the state of art technology. 

PACS numbers: 03.75.Lm,05.30. Jp,03.75.Dg,34.10.+x 



INTRODUCTION 

The high controllability of the experimental param- 
eters of ultracold atomic systems renders them versa- 
tile candidates for the implemention of quantum infor- 
mation protocols. In particular, the ability of tailoring 
traps of various geometry [l|, Q and of tuning the inter- 
atomic interactions [sj makes ultracold atomic systems a 
unique playground for applications in quantum technolo- 
gies. Recently, the Bosonic Josephson Junction (BJJ), 
a system of two-mode ultracold bosons, has received a 
large theoretical and experimental interest. 

It has been proposed 111 [l^ and experimentally demon- 
strated [8l4lO| that such a system allows for the creation 
of atomic squeezed states. These non classical states can 
be used in an interferometric protocol to improve the 
phase sensitivity reducing it below the shot-noise limit 
1/Vn [13|. as recently demonstrated in BJJs by C. 
Gross ei aZ [9|. A further enhancement of precision in 
atom-interferometry has been predicted to be reached by 
the use of macroscopic superpositions of atomic coher- 
ent states. Such states allow in principle for a scaling 
of the phase uncertainty as 1/iV, which corresponds to 
the "Heisenberg limit" [6,[l4|. Their experimental real- 
ization would be hence of primary importance both for 
applications in metrology, and for the fundamental inves- 
tigation of the boundary between classical and quantum 
dynamics [isj . 

These states however have so far been realized rather 
with a mesoscopic number of particles up to few tens in 
systems of ions |16| , nuclear spins 17 [ or photons in a cav- 
ity [18|. Several proposals to create macroscopic quan- 
tum superpositions in BJJs have been advanced ilia- 



fragility with respect to decoherence, induced by par- 
ticle losses J2 3il. p hase noise , collisions with ther- 
mal atoms [25I |26|. interaction with the electromagnetic 
field [27^, and random fluctuations of the trapping po- 
tential [28]. Therefore, it is of great interest to provide a 
protocol to create macroscopic quantum superpositions 
in the shortest time possible, i.e. before that decoher- 
ence becomes effective. The optimization of the fidelity 
with which such states can be formed in a given time, 
this latter being chosen as small as possible, represents 
another important task. Some intuitive techniques have 
been developed over the past few years to fulfill such ob- 
jectives (5l.[l9l- 21. laH^. However, these approaches are 
not able to reach the physical limits of the best possible 
performance in terms of duration and efficiency. These 
limits can be established by using optimal control the- 
ory [31h33| . Optimization techniques have already been 
applied in the context of ultracold atoms (see [sJllll to 
cite a few) to improve, e.g., the performances of an in- 
terferometer [3^. In this paper we apply such a method 
to provide a protocol for the efficient creation of macro- 
scopic superpositions in BJJs. In particular, we show 
how an almost perfect macroscopic superposition can be 
created in a relatively short time by means of a numeric 
optimal control approach. 

Before doing this, it is instructive to consider a geo- 
metric optimal control approach 
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22| . The realization of superposition states with a macro- 
scopic number of particles is challenging because of their 



This method 

provides a generalization of the proposal of Micheli et 
al reported in Ref. [H, and allows to obtain an estima- 
tion for a bound on the time of formation of macroscopic 
superpositions in the bosonic Josephson junction. In a 
summarized way, geometric optimal control is a vast do- 
main where the optimal control problems are solved by 
using tools of geometry and Hamiltonian dynamics. Sur- 
prisingly, while optimal control is known in quantum me- 
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chanics since the eighties, only few results of g eometric 
control have been published up to date iO, 41|. Recent 
advances in this domain have permitted to attack prob- 
lems of increasing difficulty. In particular, some of the 
authors have applied this method with success in fun- 
damental quantum control problems 4^ and in Nuclear 
Magnetic Resonance [i^. Due to its geometric frame- 
work, this method is however intrinsically limited to sys- 
tems with few degrees of freedom. In our example, we 
describe the BJJ in the two-mode approximation by the 
two-site Bose-Hubbard model, and we consider its semi- 
classical limit, which is valid when the number of parti- 
cles is large. In this limit, the system can be described 
in terms of two classical conjugated variables, which can 
be related to the polar and azimuthal angles of a Bloch 
sphere of radius where N is the total number of 

atoms. Hence the tools of geometric optimal control the- 
ory can be applied. Within this approximation, we show 
how to reach in minimum time a quantum superposition 
state. Due to the semi-classical approximation, the effi- 
ciency of the optimal solution is however limited in the 
original quantum domain. 

We then determine the solution of the initial quan- 
tum problem by using a purely numerical approach, the 
monotonic convergent algorithm. The efficiency and the 
flexibility of such algorithms have been demonstrated in 
a large number of studies [4J, |45| . In order to guide the 
numerical optimization, we use a coupling between the 
geometric and the numerical approaches in the sense that 
the geometric solution is used as a solution to initialize 
the numerical algorithm. The duration of the numerical 
solution necessary to produce with efficiency a superpo- 
sition state is shown to be about ten times larger than 
with the geometric approach. We analyze the robustness 
of this factor with respect to the number of particles. 

The paper is organized as follows. We introduce in 
Sec. □ the model system, i.e. the bosonic Josephson 
junction in the two-mode approximation. We also re- 
call how to construct the corresponding classical limit 
when the particle number is large, and we set the control 
problem under study. We then briefly recall how to pro- 
duce macroscopic superpositions within the protocol of 
Ref. [HI . Section □ focuses on the numerical results of the 
geometric and purely numerical methods, a comparison 
with the existing results is also made. Finally, we ana- 
lyze the robustness of the optimal solution against atom 
number fluctuations. In Sec. „ we discuss the possible ex- 
perimental implementations of our o ptim al solutions with 
the state-of-the-art technology [1, |9|.|46| . Conclusion and 
prospective views are given in a final section. Section □ of 
the appendix is devoted to the presentation of the geo- 
metric and numerical optimal control approaches. Some 
technical computations are also reported in Sec. □ and □ 
of the appendix. 



THE BOSONIC JOSEPHSON JUNCTION AND 
THE SETTING OF THE PROBLEM 

The model system 

The Bosonic Josephson Junction can be described in 
the two-mode approximation (assuming fixed the total 
number of particles N) by the spin-Hamiltonian 



H = xJl + 5J, - nj, , 



(1) 



where the angular momentum operators Jx, Jy, and 
are related to the annihilation operator a, of an atom in 



the mode j = 1,2 respectively by Jx 



J, 



'i{a[a2 



4«i)/2, 



and Jz 



alai)/2, 
a^a2)/2 



{a\a2 
n = {a\ai 

(number imbalance operator) [47l |4^. This model de- 
scribes both a system of cold atoms in two distinct hy- 
perfine states, trapped in a harmonic potential (internal 
BJJ) and a system of bosons confined in a double 

well potential (external BJJ) The interaction con- 
stant, X = Ui + U2 — 2Ui2, related to the atom-atom 
interaction energy Ui in the mode i and to the cross- 
interaction U12, can be tuned to some extent, by exploit- 
ing e.g. Feschbach resonances The cross-interaction 
term U12 is usually taken to be zero in the external BJJ, 
because of the spatial separation of the two modes. The 
parameter 6 is related to the energy imbalance between 
the two modes and to the interaction imbalance, and is 
defined by 5 = {Ei - E2) + {N - l){Ui - U2)/2. In the 
external case the Rabi coupling i7, which represents the 
tunneling term, is usually positive, i.e. £7 > 0, and can 
be tuned by changing the height of the barrier separating 
the two wells. A negative coupling Q. can be engineered 
by applying a drive to the barrier height 'H^. Such a 
coupling is however far more efficiently controllable in 
the internal setup, where by modulating amplitude and 
phase of resonant microwave and radiofrequency fields 
shone on the atoms, both amplitude and sign of VL can 
be tuned instantaneously with respect to the other time 
scales of the problem \^ (see also the concluding section 
for a discussion about the experiments). Hence, to de- 
velop our control protocol we will keep fixed the param- 
eter X with 5 — and use £7 as control field, having in 
mind the internal BJJ setup which appears more suitable 
for the experimental implementation of the protocol. 

As initial state of the dynamical proposal described in 
the following we will take an atomic coherent state where 
all of the atoms occupy the same one-particle state [sij 
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with a = tan(6'/2) exp(— i(/i), and \n) the Fock state sat- 
isfying 



Jz\n) = n\n). 



(3) 
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Each coherent state ([2]) can be represented on 
the Bloch sphere as a circle, the center of which 
has coordinates given by the expectation value of 
the angular momentum operator {d^(j)\J\9^(j)) — 
A^(sin0cos^, sin0sin0, — cos0)/2. Since the quantum 
fluctuations of the angular momentum operators in each 
direction tangential to the Bloch sphere in the point 
{9,4>\J\9,4>) are given by ^/{Kjif = yjN/A, as an or- 
der of magnitude for the radius of the circle we can take 
CT = y/N . Starting from the Heisenberg equations for J, 
one obtains by using a semi-classical limit Q that the 
angular coordinates {9, (j)) satisfy the system 



i.e. the superposition of the two coherent states located 
on the equator of the Bloch sphere along the x- axis. 
Among all the functions uj{t) allowing to reach such tar- 
gets, the optimal solution will be the one minimizing a 
given cost, such as the duration of the control (with a 
bound on the amplitude of the control field) or the en- 
ergy of the field (with a fixed control duration). In this 
framework, optimal control theory can be viewed as a 
vast machinery aiming at solving such problems either 
analytically (or with a very high numerical precision) for 
geometric methods, or numerically for the optimal con- 
trol algorithms. 
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— il sin0 
S — Nx cos 9 — Q. cot 9 cos ( 



which can be written as follows 
XN 



—uj sm <p 
A — 2 cos 9 — UJ cot 9 cos < 



(4) 



(5) 



by introducing the parameters A — 25 /{x^) ^-^d uj ~ 
2Vt/{xN). In the following, we consider the situation 
where 5 = and the scalar factor is assumed to 

be constant. The system is then controlled by only one 
parameter, namely uj. Equation ([5]) can also be written 
in a more compact form as 



F + ujG 



(6) 



where F and G are two vector fields of coordinates 



xN 



(0,-2cos^ 



and ^(- 



cot 9 cos < 



The dynamics generated by the Hamiltonian ([T|) can 
produce macroscopic sup erpositions of coherent states at 
specific times 0, S3- In order to create such super- 
positions in an optimized way, we formulate this problem 
as an optimal control problem. We first assume that the 
system is initially in the state |7r/2, tt), the goal being to 
design a "field" uj{t) (we follow here the control terminol- 
ogy) such that the system reaches a superposition state. 
We will consider as target states perfect macroscopic su- 
perpositions of coherent states ©, commonly referred to 
as "Schrodinger's cat state". For example, the state 



|Cati) = -i=(|^ = O) + |0 = ^)) 



(7) 



is the superposition of the two coherent states on the 
poles of the Bloch sphere. This state is also know as a 
"NOON cat state", because of its equivalent expression 
on the basis of the mode occupation |Cati) = -^(1-^, 0)-f 
|0, N)). The rotation of the previous state by 7r/2 around 
the y axis leads to the "phase cat state" 



I Cats) 



1 



V2 



(|^? = 7r/2,</) = O) + |0 = V2,0 = 7r)), (8) 



Creation of macroscopic superposition in the 
proposal of Micheli et al 

In Ref. [1] , a proposal for the creation of a NOON state 
has been given, based on a semi-classical argument. We 
briefly recall here this method, in the direction of which 
the experiments of Ref. [i^l are performed. 

In the classical model, the dynamics can be visual- 
ized by means of trajectories on the Bloch sphere pa- 
rameterized by the coordinates 9 and (j) (see Fig. [H and 
also Refs. Slllii). For w = 2 a bifurcation occurs in 
the model, and the fixed points of the system change: 
in the Rabi regime (w > 2) i^o = (^' = 7r/2,(/i = 0) 
and = {9 = 7r/2,(/) = tt) are the two stable fixed 
points, while in the Josephson regime (a; < 2) be- 
comes unstable and two new stable fixed points at F± = 
{cos 9 = - (w/2)2,(?!) = tt) appear In the 

latter regime, trajectories in which the number imbal- 
ance can not be reduced to zero are allowed (macroscopic 
quantum self-trapping). For each value uj < 2, there ex- 
ists hence a special eight-shaped trajectory which sepa- 
rates the macroscopic quantum self-trapping trajectories 
from the oscillations in which the number can take the 
zero value. This special trajectory, which passes through 
Ftt, is called the separatrix and can be used to produce a 
quantum macroscopic superposition state. In particular, 
ioi UJ — 1 the two poles of the Bloch sphere belong to the 
separatrix. The initial state of the quantum dynamics 
is a coherent state centered in F^. As recalled in Sec. „ 
the width associated to its fluctuations is fi = ^/N. In 
a semi-classical picture, such a state can be viewed as a 
cloud of points, a half of them being in the upper hemi- 
sphere of the Bloch sphere and a half in the bottom one, 
each of them evolving according to the classical trajecto- 
ries. As a result, for cj = 1 the initial wave packet evolves 
along the separatrix by splitting into two outgoing parts, 
one going to the north pole and the other to the south 
one. From a control point of view, the solution of Ref. Q 
consists of choosing a field of constant amplitude, i.e. a 
bang a; = 1. The time that it takes for a point initially 
in (9 = 7r/2,0 = tt) to travel along the separatrix and 
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reach one of the poles of the Bloch sphere is [5| 

XT.. (9) 

which can be taken as an estimation of the time of 
formation of the macroscopic superposition. Note that 
this time is much smaller than the time required to 
form a macroscopic superposition of phase states by the 
quenched dynamics of the BJJ (i.e. the dynamics driven 
by the interatomic interactions only) which is given by 
= 7r/2 [igj . providing therefore a speed up with re- 
spect to the use of that protocol. One may then check the 
efficiency of the semi-classical process presented above by 
performing a quantum calculation of the time evolution 
of the system under the Hamiltonian ([1} with uj ~ 1 and 
6 — 0. We compute the projection of the state at time 
Tc on the cat state given in Eq. © i as introduced in 
Appendix □ Such a calculation leads to a projection of 
Pi = 0.1394, with the parameters x = 1, ^ = 300. 

CONTROLLED CREATION OF MACROSCOPIC 
SUPERPOSITIONS 

Control of the semi-classical model 

We propose to revisit the analysis presented in the pre- 
vious section in the framework of geometric optimal con- 
trol theory. In this section, we apply this approach, i.e. 
the Pontryagin Maximum Principle (PMP) and the cor- 
responding numerical tools (see Sec. □ of the appendix 
for a presentation of this approach) on the semi-classical 
model on the sphere whose dynamics is governed by Eq. 
([5]). The control problem consists of creating a Cat state 
in minimum time with a field bounded by the value m. 
The role of this bound will be also analyzed in this para- 
graph. We denote by Tmin the minimum time required 
to reach the target state. 

We consider as initial classical state a point on the 
Bloch sphere at a distance a oi {9 = 7r/2, 4> = tt), cor- 
responding to the extremum point on the uncertainty 
circle of the coherent state 16* = TT/2,(j) = tt). Without 
loss of generality, we can choose this point in the up- 
per hemisphere. With this reduction, the semi-classical 
control problem becomes of dimension two and all the 
tools presented in Sec. □ of the appendix can be used. 
The application of the PMP leads in the spherical coor- 
dinates (9, (f) to the pseudo-Hamiltonian if, analogous 
to Eq. ^ 

H ^p- (F + ujG). 

Introducing the components {pe,Pcji) of p, the Hamilto- 
nian H reads as 

H — —2 cos 9p^ — ujpg sin (f> -— tup^ cot 6 cos 0. 



Straightforward calculations following the line presented 
in Sec. □ of the Appendix lead to det[F, G] = — 2cos6'sin0 
and det[G, [F, G]\ = 2 sin^ 0(1 - cos^ 9) - 2 cos^ 9. The 
singular set S depicted in Fig. [T]is therefore given by 

(j) = arcsin[± cot 0], (10) 

while the singular control field is equal to 

u!s = sin COS (/i. (11) 

At this point, we can compute the optimal sequence 
to reach in minimum time the north pole of the Bloch 
sphere. By symmetry of the dynamical equations, the 
point of the lower hemisphere symmetric with respect 
to S of the initial state will reach simultaneously the 
south pole. This classical simultaneous control leads in 
the quantum domain to the creation of the superposition 
state |Cati). In the optimal control theory framework, 
note that there is no restriction on the value of the bound 
m of the control field but simple solutions can only be 
obtained if the bound satisfies m > 1. We will focus on 
this case in the numerical examples. 

Our control protocol can be generalized to create a 
phase cat |Cat2) in addition to the state |Cati), which is 
not possible in the original proposal of Ref . [5[ . 

The results of the time-optimal control problem for 
reaching both states with various bounds show that no 
significant improvement in terms of speed-up nor fidelity 
is achieved, compared to the protocol of Ref. Ta- 
ble U] summarizes the numerical results, while a detailed 
analysis is reported here below. 







NOON 


PHASE 


m 


1 


2 


100 1 2 100 


Xt(xl0-3) 
P 

Fq/N^ 


25.9 
0.139 
0.636 


24.6 
0.122 
0.596 


23.6 25.5 24.6 23.6 
0.116 0.091 0.100 0.116 
0.587 0.514 0.547 0.586 



TABLE I: Numerical results of the semi-classical control pro- 
tocol for three different bounds, m = 1, 2 and 100. The 
control duration (xi), the projection (P) and the Fisher in- 
formation (Fq/N^) are given for the two cat states |Cati) 
(NOON state) and |Cat2) (Phase state) considered in this 
work. In the quantum calculation, the parameter N is taken 
to be A = 300 



We first begin the analysis by the NOON state |Cati). 
Different optimal trajectories are displayed in Fig. [T] for 
three bounds, namely m — 1, 2, 100. In the first situa- 
tion, for m = 1, we recover the solution of Ref. [11 which 
is only composed of a bang pulse, since the pole of the 
Bloch sphere can be reached by following the separatrix. 
More complicated solutions can be constructed when the 
bound m takes larger values. In this case, one can follow 
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the singular line from the initial point of the dynamics 
to the intersection point between the separatrix and the 
singular locus. To simplify the discussion, we assume 
here that this initial point belongs to the singular set. 
We then leave the singular set by using a bang pulse to 
follow the separatrix and to reach a pole as can be seen 
in Fig. [1] When m 1, bang trajectories are very close 
to the meridians of the Bloch sphere. A comparison be- 
tween the optimal solutions for m = 2 and to = 100 and 
the solution of Ref. 0] is displayed in Fig. [TJ With this 
solution, we reach the target in a time xTmm — 0.0236 
for TO = 100 while, from Eq. ©, x^c = 0.0259 in the 
same conditions for m = 1. The corresponding projec- 
tions on the state |Cati), obtained by calculating numer- 
ically the time-evolution of the quantum state under the 
field a;(i), are equal to Pi = 0.116 and Pi = 0.139. The 
values of the Fisher information are Fq/N'^ = 0.539 and 
Fq/N'^ — 0.636 respectively, which demonstrates that 
correlations among the component of the superposition 
exist, since Fq ~ > TV for = 300. The defini- 

tions of these quantities are recalled in Appendix □ The 
probability distribution of the Fock states is shown in 
Fig. [TJ As can be expected, such a distribution for the 
superposition created with our protocol mainly presents 
two peaks, close to the one of a perfect NOON state (see 
Appendix n for details). 

The same arguments can be used to describe the opti- 
mal trajectories reaching the state |Cat2). Figure [2] shows 
that the structure of the solutions are very similar for the 
different values of the bound to. The dynamics follows 
the singular line up to the separatrix, where a bang pulse 
is used to reach a point of the equator. Here, note that 
the sign of the bang pulse is different from the one of the 
first example of Fig. [TJ Numerical results comparable 
with the state |Cati) are obtained in this case, as can be 
seen in Table [TJ The numerical results for the projection 
on the target state and the Fisher information as a func- 
tion of time are reported in Fig. [5J Fringes appear in the 
Fock state distribution of the final state as explained in 
Appendix □ 

The time Tmin can be estimated analytically for any 
value of TO by integrating the classical equations along 
the singular line. For m — > +cx), we have 

2 dx 

Jm ^sin^ x{l + sin^ x) 

where 4'{Q) is the coordinate of the initial point on the 
singular set. Hence Tmin is inversely proportional to the 
total number of particles as in the case of solution of 
Ref. 0, from which Eq. ([T2j) differs by a numerical fac- 
tor. The computation leading to Eq. ([T2t . as well as a 
numerical comparison between Tc and Tmin , are detailed 
in Appendix n 

To conclude, we stress that in these different time- 
optimal-computations both the minimum time of forma- 
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FIG. 1: (Color online) (left) Plot in the (0, (^) plane of the 
optimal trajectory in the semi-classical model for the bounds 
m — 1, m = 2 and m — 100 in blue (dark), green (light gray) 
and red (dark gray) respectively. The target state is the cat 
state |Cati). The dashed line indicates the position of the 
singular set. The dashed blue, red and green lines represent 
the position of the separatrix for the different bounds, (top 
right) Evolution of the corresponding control fields a; as a 
function of the dimensionless time x^- (bottom right) Fock 
states distribution P{n) = |(n|^/)(f))p as a function of n for 
t — 0, t = Tmin and the target state |Cati). In the quantum 
calculation, the parameter A'^ is taken to be A'^ = 300. 



tion of macroscopic superpositions and the respective fi- 
delities are comparable to the ones obtained with the 
method of Ref. [5| (the numerical comparisons are sum- 
marized in Table [I|) . This shows the efficiency of the so- 
lution of Micheli et al for controlling the semi-classical 
model, and suggests that a time of order of Tc is the phys- 
ical bound on the time of formation of a macroscopic su- 
perposition in a BJJ. Note that this non trivial result has 
been derived from the tools of geometric optimal control 
theory. To get better results in terms of fidelity, another 
approach able to tackle the fully quantum character of 
the problem has to be used, which we will address in the 
following section. 
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FIG. 2: (Color online) Same as Fig. [T]but for the target state 
ICata). 

Control of the quantum model 

We apply in this paragraph the monotonic algorithm 
presented in Sec. □ of the Appendix to maximize the 
projection onto the target state at time t — T . We will 
consider different control durations, namely T = Tc, T — 
bxTc and T = 10 x Tc , which are multiple of the minimum 
time Tc for a bound m — 1. These different tests will 
allow us to determine the time required to produce a cat 
state with efficiency. 

The Hamiltonian ([l} of the system is written as 

H = Ho + ljHi, (13) 

where Hq = xJzi ^1 = ^^^x ^^'^ ^ now the new 
control parameter. In the following computations, we 
choose as parameters of the quantum system x = 1 and 
N — 300. The parameters of the algorithm are respec- 
tively taken to be (A = lO"*', 771 = 772 = 10^) for T = 
and (A = 5 X 10~^, ?7i = 772 = 2) for the two other cases. 



The different computations are presented in Fig. [3] and 
El for two durations T = 5 x and T = 10 x T^. We 
consider that the initial field of the algorithm is the so- 
lution of Ref. 0, i.e. a constant field u}(t) = 1 in the 
interval [0, T]. Due to the structure of the control prob- 
lem, note that the algorithm does not converge towards 
an efficient solution for any initial condition. Taking as 
initial field the solution uj(t) corresponding to higher val- 
ues of the bound in the geometrical protocol would not 
change significatively the result because of the proximity 
of this solution with the one of Ref. Q . 2000 iterations 
are used in the first example, 1000 in the two others. 

We obtain that for the case T = 5 x and T = 10 x 
the target macroscopic superposition can be created with 
a very high fidelity, respectively larger than 0.88 and 0.98, 
corresponding to Fq/N"^ = 0.951 and Fq/N'^ = 0.997. 
Note that for T — Tc, a projection of 0.2548 is reached, 
which is better than the result provided in Ref. 0], but 
at the price of a more complicated solution. This result 
was expected in the sense that the solution by Micheli et 
at has been derived in the semi-classical limit without 
taking account the quantum character of the problem. 
A visualization on the Bloch sphere of the state created 
is provided in the second panel of Fig. The same 
computation has been done for the target state |Cat2) 
for T = 10 X Tc and the parameters A = 5 x 10~^ and 
771 = 772 = 2. As initial condition for the algorithm, we 
consider the bang-bang optimal solution with a total du- 
ration increased by a factor of 10 (we keep constant the 
relative time of the two bangs). We recall that a bang 
pulse is a field of maximum intensity and a bang-bang 
sequence, the concatenation of two bang pulses of differ- 
ent signs. In Fig. El note the bang-bang structure of the 
final solution constructed by the algorithm. This exam- 
ple shows that the algorithm has been guided towards a 
particular mechanism by the semi-classical computation. 
The different numerical results are listed in Table El A 
visualization on the Bloch sphere of the state created is 
provided in the third panel of Fig. |4l 



NOON PHASE 

n 1 5 10 1 5 10 
P 0.255 0.880 0.994 0.245 0.903 0.989 
Fq/N'^ 0.650 0.951 0.997 0.632 0.946 0.996 

TABLE II: Same as Table U but for the quantum protocol. 
The parameter n represents the ratio of the control duration 
over the time Tc (see the text). 



In these control problems, we observe that a minimum 
time required to produce with efficiency a cat state (i.e. 
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FIG. 3: (Color online) Plot of the evolution of the projection 
of the state IV'(i)) onto the target state |Cati) (middle), of 
the control field lo (top) and of the Fisher information Fq 
(bottom). The control duration is 5 x T^. In the top panel, 
the horizontal solid line is the solution of Ref. Q, which is 
taken as a trial field of the algorithm. 




FIG. 4: (Color online) Plot of the projections on the Bloch 
sphere, i.e. the Husimi function Q{9,cj>) = of the 

initial state |7r/2, tt) (left), the final states with the geometric 
solution where m = 1 and T — Tc (middle) and with the 
monotonic algorithm (right) where T — 10 x Tc- 




2 4 6 8 10 



t/T 

c 

FIG. 5: (Color online) Same as Fig. [3] but for T = 10 x T, 



FIG. 6: (Color online) Same as Fig. [5]but for the target state 
I Cats). 

with a projection larger than 0.99) is of the order of 
10 X Tc- For the NOON state, we have analyzed the 
variation of this scaling factor with respect to the parti- 
cle numbers TV. Numerical simulations revealed that for 
N ranging from 50 to 700 the target state is created with 
a very high precision for this duration, as can be seen 
in Table IIIII (to simplify the discussion, we only indicate 
the projection onto the target state and the Fisher in- 
formation for a time of 10 x Tc, even if we have tested 
other control durations). This shows that 10 x Tc is a 
relevant time for this control problem, at least in the ex- 
perimentally range of particle numbers considered. For 
low values of A'^, it roughly corresponds to the minimum 
time where an efficient control field can be designed. For 
all the values of N considered, the optimal solution is 
very close to a constant solution as shown in Fig. [5] with 
N = 300. Note that, when N > 500, we have numerically 
observed that lower times of the order of 7 x Tc can be 
chosen, but at the price of a more complicated structure 
for the field. This seems to be a general property of these 
quantum computations, i.e. the factor multiplying Tc to 
consider to produce a superposition state decreases as N 
increases. The determination of the exact dependence of 
this factor with the particle number is a hard task, since 
by construction the numerical algorithm does not allow 
to seek for the optimal time. 



In a final step, we have analyzed the robustness of the 
optimal control field as a function of A'', i.e. against the 
atom number fluctuations. We choose as reference con- 
trol field the one corresponding to iV = 700. We have 
determined for TV ranging from 650 to 750 the projec- 
tion and the Fisher information obtained at = 10 x Tc 
with this control field. This point is illustrated for the 
NOON state in Fig. [71 We observe that the efficiency of 
the process quickly decreases when the particle number 
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N P F/N^ = 10 X TcjN) 

50 0.9995 0.9997 1.1976 

100 0.9995 0.9996 0.6634 

200 0.9997 0.9998 0.3688 

300 0.9992 0.9997 0.2594 

500 0.9998 0.9999 0.1659 

700 0.9996 0.9999 0.1233 

TABLE III: Numerical results of the optimal control compu- 
tations for different atom number A^. The control duration 
{xt = 10 X Tc), the projection (P) and the Fisher information 
(Fq/N^) are given for different values of A^. 




FIG. 7; (Color online) (top) Plot of the projection P on the 
target state as a function of the variation of atom number 
A N for the optimal control field computed at A = 700 (see 
the text). The control duration is taken as 10 x Tc. (bottom) 
Same results but for the Fisher information. 



changes. For instance, the average projection is already 
of the order of 0.5 for a difference of 10 particles. This 
behavior is not surprising in optimal control where the 
optimal solutions are known to be very efficient (even for 
very complicated problems) , but are known to lack of ro- 
bustness. A method based on the simultaneous control 
of different systems can be used to improve this latter 
property (see, e.g., Ref. [E^l)- This approach, which can 
be very time consuming in terms of numerical compu- 
tations, requires in addition to consider longer optimal 
solutions than the ones used in this paper. This means 
that the limit time of 10 x Tc is no more valid for de- 
signing an optimal solution robust with respect to large 
variations of N. Due to the complexity of this method, 
such a computation goes beyond the scope of this article 
and has not been undertaken. 



EXPERIMENTAL IMPLEMENTATION OF THE 
PROTOCOL 



We now analyze to which extend the fields correspond- 
ing to the solutions found could be implemented in real 
experiments with the present state of the art technology. 
Let us discuss first the amplitude of the control. ^From 
Fig. [5J we see that the solution for the field uj, leading to 
the formation of a NOON state with a very high fidelity in 
a time 10 x Tc, is roughly comprised between —1 and 1.5. 
Having in mind the internal BJJ setup (and in particu- 
lar the experiments of Ref. (ol. |46|). typical experimental 
bounds on the parameter J7 are < < 27r • 2 kHz, and 
a typical value for x is x ~ 27r • 0.13Hz. Fixing such a 
value for the interactions translates into < < 102.6, 
which largely comprises the values required by our so- 
lution. Incidentally, we remark that even the geometric 
solutions presented in Section □ would be implement able, 
the largest value of w considered there being w — 100. 
Furthermore, the control field can be switched fast com- 
pared with the other time scales of the experiments. Note 
that it would be also possible to include in our protocol 
some spectral constraints on the control field [45*1. Hence, 
in ideal conditions it would be possible to implement our 
control protocol. 

However, in realistic conditions the experiments are af- 
fected by the presence of noise, which induces dissipation 
and decoherence. The main sources of noise in the BJJ 
system are represented by particle losses Q and phase 
This latter is due to stochastic fluctuations of 



noise 



the energies of the two modes of the BJJ. Other sources 
of noise are discussed in Refs. 25l427l|. A quantitative 



discussion of the effect of decoherence on the creation of 
superposition states would require to treat on the same 
footing in our control protocol the unitary dynamics lead- 
ing to the formation of the cat state and the decoherence 
sources. This point goes beyond the scope of this work, 
but represents an interesting perspective, which may be 
useful for the experimental quest of a macroscopic super- 
position. For a qualitative discussion, we can take the 
order of magnitudes of the decoherence rates of a process 
in which the macroscopic superposition is taken as initial 
state. We assume that this state evolves under the effect 
of noise only. In this case, we may distinguish the dissipa- 
tion time Tdiss, which characterizes the relaxation of the 
components of the superposition, from the decoherence 
time Tdoc, i.e. the time at which the coherences of the 
superposition are washed out [s^ . Issf . In general these 
two time scales are well separated. An exception is given 
by the decoherence induced on phase cat states by phase 
noise, for which the decoherence rate does not depend 
on the total particle number [24] (this is not the case for 
the NOON cat, for which t^^^'' oc N under phase noise). 
Losses processes can be one-body (due to scattering with 
impurities), two-body (spin- relaxation losses, due to col- 
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lisions of two atoms which change their spin state, free- 
ing a high kinetic energy and ejecting them out of the 
trap), or even three-body (when three atoms cohide, two 
of them form a molecule, ejecting the third out of the trap 
(60|). The decoherence times are respectively t^~^^ oc N 



for one body- losses, 



(-1) 

dec 



A'^^ for two-body losses and 
"^doc^^ oc N'^ for three body losses (see e.g. Ref.fHot for 
the three body case). Hence, if one-body losses and/or 
phase noise are affecting the system, since the time of for- 
mation of the cat states within the control protocol is of 
the order of 1/7V and the decoherence rate grows linearly 
with the atom number, it is not obvious to determine the 
most convenient atom number. An analysis close to this 
spirit in the context of the creation of squeezed states in 
the presence of particle losses has been carried out by Li 
et al in Ref. [23|. If two and three-body losses are non 
negligible, since the decoherence rate grows respectively 
with N'^ and the experiments should be limited in 
practice to small atom numbers. 

Let us consider a gain the order of magnitudes of the ex- 
periment in Ref. (jl l46l| , where the initial particle number 
is N Ki 400 (we have used N — 300 in our simulations). 
The loss rate depends on the interaction strength, and is 
typically of Tdiss ~ 0.1s for x ~ 27r x 0.13 Hz [ig. This 
corresponds to loosing about 15% of the atoms in about 
20 ms. Let us suppose that one-body losses are the dom- 
inant noise process. The time for the efficient generation 
of a NOON state is 10 x T^, resulting in 10 x xTc = 0.259, 
i.e. Tc — 0.317 s, which is already of the same order of the 
dissipation time Tdiss- Furthermore, the decoherence rate 
would be in this case Tdcc ^ Tdiss/SOO ^ 0.33ms, which is 
a thousand time smaller than the time 10 x Tc- There- 
fore, in order for our protocol to be implementable with 
the experimental parameters \^ , the dissipation rate 
should be lowered by a factor 1000. In this framework, 
an analysis of the optimum atom number to use in the 
presence of decoherence could yield useful indications. 



CONCLUSIONS 

In this paper we have designed different control fields 
by using geometric and numerical optimal control tech- 
niques in order to create quantum superpositions in a 
bosonic Josephson junction. We have generalized the 
solution proposed by Micheli et al. in Ref. [HI for the 
semi-classical model. The minimum time Tmin of for- 
mation of macroscopic superpositions via our geometric 
control protocol has been shown to be of the order of 
1/A^, improving the result of by a numerical constant 
(see Appendix □ for the detailed calculation) , with a fi- 
delity of the same order. This also shows the efficiency 
of the solution by Micheli et al., which is very close to the 
time-optimal solution in the semi-classical framework. 

An almost perfect fidelity can be reached in longer 



times of the order of 10 x Tmim with a more complicated 
field solution. We have checked that this duration is suf- 
ficient for creating very good superpositions for a large 
range of particle number N , the fidelity further improving 
at increasing N . Note that such a duration can still lead 
to a significative speed-up with respect to the protocol 
based on the quenched dynamics of the BJJ [l9|, where 
the characteristic time of formation of a two-component 
superposition of phase states is 7r/(2x), independent on 

N m- 

Decoherence effects may change significatively the 
state crated with the designed control field, and a careful 
analysis of the control sequence to be used in the presence 
of noise stems as an important extension of the present 
work. 



Optimal control methods 
Geometric approach 

We present in this section the geometric optimal the- 
ory which allows us to solve the optimal control problem 
in the semi-classical limit. In this work, we consider a 
system belonging to a general class of optimal control 
problems on a two-dimensional manifold (here the Bloch 
sphere) with a single control field. Powerful mathemati- 
cal tools that we briefly describe below have been devel- 
oped in this case [sjj . 

We consider a system governed by a differential equa- 
tion of the form 

f = F{x) + uG{x), 

where a?' is a two-dimensional state vector, F and G two 
vector fields and u a control field satisfying the constraint 
|u| < m, where m is a bound on the control. Starting 
from a point xq, the goal of the control is to reach in 
minimum time a target Xf. We assume that a solution 
exists, i.e. that the target belongs to the accessibility 
set of the initial point xq. We solve this time-optimal 
control problem by applying the Pontryagin Maximum 
Principle (PMP), which is formulated from the pseudo- 
Hamiltonian 



H = p- (F + uG), 



(14) 



where p is the adjoint state. The Hamiltonian trajecto- 
ries which are solutions of 

. dH ^ . dH 

X = -g^[X,P,V), P = --Q^(X,P,V), (15) 

where the control field v is determined from the maxi- 
mization condition 



H{x,p,v) 



max H{x,p, u), 

\u\<m 



(16) 
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with H > 0, are candidates to be optimal. Such trajecto- 
ries are cahed extremals in the control community. Other 
geometric tools have to be used (see [s^ for a recent 
overview on this point) to select among the extremals, 
the ones which are effectively optimal. 

The maximization condition (jl6p can be solved from 
the switching function $ given by 



(17) 



Using the maximization equation, one deduces that an 
extremal curve {x{t),p{t)) is composed of a sequence of 
arcs 7+ , 7_ and 7s . Regular or bang arcs 7± are obtained 
when sign[$] = ±1. In this case, the control field is given 
by u = TO X sign[$]. When the function $ takes a zero 
value and when this zero is isolated, the field switches 
from ±m to =Fto. We encounter a singular situation and 
singular arcs if $ vanishes on an interval [to,ti]. In this 
case, the control field cannot be directly determined from 
the maximization condition, but from the fact that $ — 
$ = •••= on [to, ii]. A straightforward computation 
leads to 

where [G, F] is the commutator of two vector fields de- 
fined by 

[G,F]=WG-F-G- WF, 

VG being the Jacobian matrix of the vector G. One then 
deduces that the system $ = $ = admits a non-trivial 
p solution on the singular set S defined by 



{f ; det(G, [G, i^])(x) = 0}, 



(18) 



where the vector fields G and [G, F] are parallel. The 
singular arcs 7^ are located in the set S. The singular 
control field Us can be computed from the second deriva- 
tive of $ 

$ = p • [G, [G, F]] + u,p- [F, [G, F\] = 0, 
which leads to 

P-[G, [G,F]] 



P-[F, [G,F]] 



(19) 



Note also that the optimal solution can follow the singu- 
lar lines only if the control field is admissible, that is if 

\Us{x)\ < TO. 



Monotonically convergent algorithms 

Monotonically convergent algorithms are a standard 
approach to solve the optimality equations in quantum 
mechanics. Such algorithms have been described in detail 



elsewhere 4J, |45|. Here, we only summarize the main 
theoretical aspects to apply these methods (iHj . 

We respectively denote by \(j)o) and \(j)f) the initial and 
target states of the control problem and we assume that 
the Hamiltonian of the system can be written as 



H = Ha + u{t)Hi, 



(20) 



where u{t) is the control field. We consider the following 
cost functional which penalizes the energy of the control 
field 



J=|(0/|V'(T))|2-A / uit^dt 



(21) 



where A is a positive parameter. The value of A expresses 
the relative weight between the projection onto the target 
state and the energy of the field. The control duration T 
is assumed to be fixed here. In order to satisfy the con- 
straint that the state \tp{t)) is solution of the Schrodinger 
equation, we introduce the augmented cost J 



j= |(<^/|^(r))|^-A / uitydt 

Jo 



where |x(t)) is the adjoint state and Im denotes the imag- 
inary part of a complex number. Using the fact that the 
variations of J with respect to \ip(t)), \xit)) and u are 
equal to 0, one arrives at a differential system of the 
form 

(z|-7?(t))|x(t))-0, 

with the initial conditions |?A(0)) = |0o) and |x(^)) = 
10/). The control field u is given by 

u{t) = -2 X ImmT)\cf,f){x{t)\H,m)]. (22) 

This set of coupled equations is solved by using mono- 
tonic convergent algorithms but other approaches such as 
gradient algorithms could be used [Hi!]. We assume here 
that the iteration is initiated from a field UQ(t). At step 
k of the algorithm, the system is described by a triplet 
(iV'fc), \Xk-i),Uk) corresponding to the cost Jk 



Jk = |(0/|^fc(r))|2- A / Uk{t)^ (23) 
Between two iterations, the variation A J is given by 

AJ^Jk+i-Jk = \{^k+i{T)\cf>f)\^-\{MTMf)\^ 
-A / {uk+i{tf -Ukit^). 
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Introducing Ak^k+i 
one arrives at 



-2Im[(V'fc+i|Xfc)(Xfc|ffi|V'fc+i)], 



AJ 



dt [{uk - Uk+i)AkM+i + (uk - Uk)AkM + 



Mul-ul+ul+i -ul)] , 



which can be written as 



AJ= / dt{Pi{t)+P2it)), 



(24) 



where 



Pi = -(wfc+i - Mfe)[A(wfc+i + Uk) + Ak,k+i], 
P2 = -{uk - Uk)[Xiuk + Uk) + Ak,k]- 

To ensure the monotonic behavior of the algorithm, we 
choose the fields Ufe+i and Uk such as the integrands Pi 
and P2 are positive. More precisely, we first determine Uk 
from Uk such that P2 is positive and then we determine 
Uk+i from Uk such that Pi is positive. This means that 
we define Uk from P2 as 



Uk-Uk = -?72[A(ufe + ilk) + Ak,k\ 



(25) 



where 772 is a positive constant. The field iik can finally 
be expressed as 



1 - Ar72 ri2 

Uk = ^ — Uk ^ — 

1 + Xri2 1 + A772 



-A 



k.k- 



(26) 



The same work can be done for the field Uk+i and Pi. In 
this case, we get 

Uk+i - Uk = -m[Huk+i + Uk) + Ak^k+i], (27) 
771 being also a positive constant, and 



1 - A7/1 , 
Uk+l = — — ; — Uk 

1 + \t]i 



■AkM+i- 



(28) 



Note that this formulation is equivalent to the one intro- 
duced in [45]. 

The structure of the algorithm can be summarized as 
follows. 

1. Backward propagation of the adjoint state \xk{t)) 
with the field Uk and the initial condition \xk{T)) — 
|(/)/). The field Uk is determined at each time 
through the formula 



2. Forward propagation of the state | V'fe+i {t)) with ini- 
tial condition \iljk+i{t)) = \4'o) and the field Uk+i 
which is computed at the same time by the formula 

(El. 



Using this procedure, it is then straightforward to show 
that the algorithm remains monotonic, that is Jk+i > Jk- 
Numerically, the algorithm is stopped when its conver- 
gence is better than a given threshold e, i.e. | Jk+i —Jk\l^ 
e. 



Measure of the coherence of the macroscopic 
superposition 

After deriving the control field w(t), we may quan- 
tify the coherence of the macroscopic superposition state 
which has been created by means of the control protocol. 
In this section, we present some quantities which will be 
used to this scope. Note that, in the case of geomet- 
ric control, such methods require the computation of the 
quantum evolution under the field uj{t). 



Fidelity 

The first method consists in computing the projection 
(or fidelity) of the final state onto the cat states intro- 
duced in Eq. 



Pi,2 = |(Cati,2|^(t))|' 



(29) 



to be evaluated at the time of formation of the superpo- 
sition. However, such a quantity does not describe the 
character of the state, and in particular its shape and its 
coherence. 



Angular momentum eigenvalues probability 
distributions 



In order to represent the macroscopic superpositions 
created with our protocol, we can use the probability 
distribution of the eigenstates of angular momentum op- 



erators in various directions |56l . |57| . These probability 
distributions carry different informations when consid- 
ering different directions and states. We consider as a 
first example the ideal NOON state given in Eq. In 
this case, the probability distribution of the eigenvalues 
associated to the Fock states, defined as the projection 
P(n) = |(n|Cati)p where \n) satisfies Eq. ([3|), is given 

by 



(30) 



This distribution displays two peaks corresponding to 
the two components of the superposition at n = ±N/2 
(see the bottom panel of Fig. [T]). The analogous Jx- 
eigenvalues distribution for the same state, i.e. Px{nx) = 
|(n2;|Cati)p where \nx) satisfies Jx\nx) = nx\nx), cor- 
responds to the profile of the NOON state when pro- 
jected on the X axis. A straightforward calculation using 
Eqs. (HI) and ^ leads to 



Px{nx) 



^(4J if 



n IS even 







if n is odd. 



(31) 



This corresponds to a binomial envelop centered in n = 0, 
with interference fringes having an unit spacing (see the 
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bottom panel of Fig. [2|) [H^ 57 1. Note that the same 
distribution for an incoherent mixture of the same two 
coherent states, such as pi = ^{\9 ~ Q,(p = Q){9 = Q,4> = 
0| + = 7r,0 = 0)(6' = TT, = 0|), would display no 
fringes. 

Since the phase cat state ([H]) is the rotation of the 
NOON state lying along the x axis, the preceding argu- 
ments have to be reversed for such a state. In particular, 
its probability distribution Pxirix) displays two peaks at 
Ux — ±iV/2, while the Fock states probability distribu- 
tion P(n) can be useful to visualize the coherences. 

Quantum Fisher information 

Additional insights about the coherence of the super- 
position state, i.e about its off-diagonal correlations, can 
be gained from the computation of the optimum quan- 
tum Fisher information i^Q = 4x maxf;(AJfi)^ [GlIHs.lHsj. 
Such a parameter was historically introduced to quai itify 
the usefulness of a quantum state for interferometry [58j , 
and it has been demonstrated that the condition Fq > N 
provides a sufficient condition for the state to be entan- 
gled [61 . In particular, for the macroscopic superposition 
states |Cati^2) given in Eqs. ([7l8l) Fq = N^, while Fq 
scales linearly with the particle number for an incoher- 
ent mixture of coherent states, such as the mixture pi 
introduced in Sec. „ In this work, this parameter will be 
used as an indicator of the presence of off-diagonal corre- 
lations in the state produced within our control protocol. 



Analytic expressions of the control durations 

We derive in this section the expressions of the du- 
rations both for the solution of Ref. [5] and the time- 
optimal trajectory. By definition of the classical control 
problem from the quantum one, we assume that the dis- 
tance between the point of coordinates O — {tt/2, tt) and 
the initial point / of the dynamics is u = \/N, this point 
belonging to the separatrix for the solution of Ref. 
and to the singular set for the optimal one. There ex- 
ist two different natural ways to define this distance, one 
corresponding to the length of the segment 01, i.e. the 
euclidian distance, the other to the length of the geodesic 
on the sphere from O to /. Denoting by {9, 0) the coor- 
dinates of /, one gets in the first case that 



-^= = \/2(l — cos (As 



the bang arc can be neglected. In this case, T^in is the 
time taken to travel from / to the point of coordinates 
{9 = 37r/4, = n/2). The dynamics is described by the 
following set of equations. 



XN 
2 ^ 



2 cos 9 — uj cot 9 cos ( 



uj = sm 

,,2 



COS 



COS ( 
2 



= sm 



ism 



where the second and third equations correspond respec- 
tively to the definitions of the singular control field and 
of the singular set. Straightforward computations lead to 

Tmin — / ; ; (32) 



2 



sin x(l + sin x) 



where </) is determined from the expression of cr. A similar 
argument can be used to compute Tc for the solution of 
Ref. along the separatrix. In this case, the dynamics 
is governed by the equations 



xN 



(—2 cos 9 — UJ cot 9 cos ( 



cos = sm f 



where coscj) = sin0 is the equation of the separatrix [Bj. 
This gives 



2dx 



1^ xNsm.x 
which can be simplified into 
2 



Tr 



xN 



\n{V2N - \/2N - 1) 



for the euclidian distance and into 



r^ = _Ain(i_^; 

xN Vl-coscr 



(33) 



(34) 



(35) 



while in the second case 



for the geodesic on the sphere. In the limit N — -|-oo, 
we recover in the two cases the expression ([9|) given in 
0. The difference between the two control durations is 
illustrated in Fig. [71 One observes that this difference 
decreases as the number of particles N becomes larger, 
which shows the efficiency of the Micheli solution in 
terms of optimality for the creation of quantum super- 
positions. 
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sin^( — ) = sin^ 6'sin^ d> + cos^ 9. 

We assume in the computation of T„iin that there is no 
bound on the control field, so that the time of travel along 
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FIG. 8: Evolution of the control duration in the semi-classical 
model (see the text) as a function of the number of particles 
A'^. The red (dark gray) and green (light gray) solid lines 
represent respectively the exact and the approximated time 
for the Micheli solution. The dashed line corresponds to the 
optimal solution with m = 1 x lO''. The small insert is a zoom 
near A*' = of the plot. 
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